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The  results  of  this  study  suggest  that  there  is  always  microscopic  chaos  present  in 
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is  different  from  the  macroscopic  chaos  which  may  be  associated  with  the  bulk 
potential  and  persists  even  for  very  large  N.  Therefore  the  Lyapunov  exponents  of 
orbits  evolving  in  TV-body  environments  do  not  decrease  as  N increases.  However 


xii 


pomtwise  comparisons  associated  with  the  power  spectra  of  the  orbits  as  well  as 
their  statistical  properties  suggest  that  the  continuum  limit  does  make  sense. 


CHAPTER  1 
INTRODUCTION 


Galaxies,  stellar  clusters,  and  galaxy  clusters  are  a few  examples  of  real 
A-body  systems  one  can  meet  in  the  context  of  astrophysics.  Of  course  a huge 
spectrum  of  A-body  systems  appears  in  every  other  field  of  physics.  It  is  desirable 
to  understand  the  dynamics  of  the  motion  of  bodies  in  such  systems,  a piece  of 
information  that  is  critical  to  our  understanding  of  the  system  itself.  Especially 
in  the  context  of  astrophysics,  real  experiments  that  involve  A- body  systems  are 
unfortunately  unrealistic  (galaxies  are  too  big  to  fit  in  the  lab!).  One  alternative  in 
our  effort  toward  understanding  is  numerical  experiments  using  digital  computers. 
Computational  physics  is  without  question  a new  paradigm  in  science,  a new  way 
to  get  answers,  wherever  analytics  or  experiments  meet  problems  too  difficult  to 
solve. 

Although  one  understands  the  interparticle  dynamical  law,  and  may  have 
some  information  about  the  mass  density  of  an  A-body  system,  it  is  usually  too 
expensive  (in  time  and  resources)  to  simulate  the  A-body  systems  as  realistically  as 
needed.  Several  codes  have  been  developed  that  can  perform  A-body  simulations. 
However  serious  resources  of  computer  time  and  hardware  are  necessary  in  most 
cases. 

Fortunately  there  may  be  an  alternative,  or  at  least  this  is  what  is  claimed. 

The  idea  is  rather  simple  in  first  view,  although  in  the  physics  terrain  there  are 
usually  subtle  dynamical  implications  associated  with  any  simplification.  Every 
system  of  A bodies  that  interacts  via  a gravitational  force  (long-range  force) 
generates  a bulk  potential  that  characterizes  the  system  as  a whole.  The  idea  is 
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easy  to  understand  if  we  think  that  the  system  of  N bodies  is  characterized  by  a 
number  density.  This  number  density  can  be  smoothed  out  to  a continuous  density. 
Then  using  Poisson’s  equation  one  is  an  integration  away  from  analytically  defining 
the  potential  (although  the  integral  cannot  always  be  solved,  one  can  still  find  the 
value  of  the  potential  numerically  for  every  point  of  interest). 

Let’s  translate  the  previous  paragraph  into  mathematics.  To  make  things 
simple  (but  without  losing  generality)  let’s  assume  that  the  total  mass  of  the 
system  is  M = 1 and  that  this  system  consists  of  N point  masses  each  of  equal 
mass  m.  Therefore  m = M/N  = l/N.  The  density  of  this  system  is  discrete  (or 
singular)  and  can  be  represented  as  a sum  of  N delta  functions,  each  of  which 
represents  one  singularity  (that  is  one  particle).  Therefore: 

p(r)  = -^'52s(r~ri)  (L1) 

i=i 

Now  this  mass  density  will  generate  a potential  the  form  of  which  can  be  found  via 
the  well-known  Poisson’s  equation: 

V2<f>  = 4t xGp  (1.2) 


where  G is  the  gravitational  constant. 

The  potential  generated  by  the  mass  density  p is  the  following: 

N 1 
7 —7 

r - r,: 


$n(  r)  = 


GM 


N 


(1.3) 


i=l 


There  is  no  approximation  until  this  point  except  that  all  the  particles  have  the 
same  mass  m.  The  equations  are  exact  and  cumbersome.  If  somebody  decides  to 
use  these  equations  to  study  the  evolution  of  the  system  he  is  doomed  to  use  one  of 
the  computationally  expensive  iV-body  codes. 


However  all  these  masses  will  generate  a potential  that  will  characterize  the 
space.  If  one  can  express  this  potential  through  an  analytical  formula  then  one  has 
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a dynamical  model  that,  in  principle,  closely  (how  closely  will  be  discussed  soon) 
describes  the  system  of  N bodies.  If  the  potential  cannot  be  expressed  analytically, 
still  one  can  express  the  density  analytically,  and  then  compute  numerically  (via 
Poisson’s  equation)  the  potential  for  any  possible  value. 

This  is  a method  for  simulating  a real  system,  but  it  has  the  complication 
that  the  iV-body  potential  cannot  always  be  translated  into  an  analytical  smooth 
potential.  However,  we  can  follow  a different  way:  There  are  many  cases  (common 
practice  in  science)  in  which  people  use  so-called  toy  models,  (models  that  are 
described  analytically  through  a simple  potential  formula  and  that  are  believed  to 
be  generic).  Both  of  these  properties  are  critical.  Simulations  of  motion  in  simple 
analytical  potentials  are  fast  and  easy  to  perform.  Generic  potentials  are  necessary 
because  one  tries  to  understand  the  physics  characterizing  more  realistic  systems 
by  studying  the  toy  model.  If  the  toy  model  is  not  generic  one  should  not  expect 
its  results  to  apply  in  general  to  realistic  systems.  Nongeneric  toy  models  are 
mathematical  ingenuities  that  have  no  application  to  nature  (although  sometimes 
applications  are  discovered  years  later.) 

Now  an  interesting  point:  Let’s  assume  one  has  chosen  an  analytic  potential 
and  has  computed  its  (analytic)  density  via  Poisson’s  equation.  He  is  now  in  a 
position  to  transform  his  model  to  an  iV-body  system.  He  is  thus  in  a position  to 
play  the  game  exactly  the  opposite  way:  from  his  analytic  model  to  go  to  an  N- 
body  system.  This  can  be  done  easily  by  using  a numerical  program  that  samples 
the  analytic  density.  What  he  produces  this  way  is  nothing  but  a representation 
of  the  analytical  model  by  a system  of  N bodies.  In  this  thesis  I deal  only  with 
time-independent  systems.  Thus  the  iV-body  analogue  of  an  analytic  potential 
will  always  consist  of  bodies  that  do  not  move  (they  have  zero  velocity).  It  is 
reasonable  then  to  “baptize”  such  a system  a “frozen- N”  system.  So  the  ultimate 
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question  is:  what  is  the  difference  between  the  dynamics  associated  with  the 
analytic  potential  and  the  dynamics  associated  with  its  frozen-N  analogue? 

Before  I analyze  the  question  even  more,  please  note  an  important  point.  It 
is  obvious  that  the  larger  A is,  the  better  the  analytic  potential  is  represented  by 
its  frozen-A  analogue.  As  A increases,  the  masses  become  smaller  and  the  typical 
interparticle  distance  also  becomes  smaller.  Thus  someone  might  visualize  that  for 
A — » oo  the  frozen- A'  potential  will  exactly  coincide  with  its  analytical  analogue. 
The  frozen-A  environment  will  be  so  dense  that  it  will  enter  the  continuous  regime, 
which  is  the  smooth  analytical  density.  This  smooth  analytical  density  is  exactly 
what  one  calls  “the  continuum  limit.” 

It  is  important  here  to  understand  the  main  difference  between  a smooth  po- 
tential and  its  frozen-A  analogue.  An  A-body  environment  is  always  characterized 
by  discreteness  effects.  The  environment  is  “grainy,”  far  from  the  smoothness  of  the 
analytical  potential.  A body  moving  in  the  frozen-A  environment  will  feel  a bulk 
smooth  potential  generated  because  gravitation  is  a long-range  force.  It  will  also 
feel  “kicks”  because  of  close  encounters  with  the  A-body  frozen  particles.  Eventu- 
ally someone  may  think  of  the  whole  motion  as  a stochastic  process  generated  by 
some  differential  equation,  where  external  potential  and  a random  force  play  their 
different  roles,  and  collectively  determine  the  dynamical  evolution.  (We  actually 
tried  that!) 

So  the  main  question  is  how  the  dynamics  in  a smooth  potential  and  in  its 
frozen-A  analogue  compare.  Before  I analyze  this  question  further  it  is  imperative 
to  describe  the  existing  information  relative  to  this  question. 

First  a numerical  result:  It  has  been  known  since  at  least  1964  that  motion 
in  an  A-body  system  is  chaotic  [1],  It  is  very  important  here  (and  it  will  become 
quite  clear  later  why)  to  emphasize  that  by  chaotic  I mean  sensitive  dependence 
on  initial  conditions  (the  classical  definition  of  chaos).  It  is  chaos  in  this  sense 
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that  can  be  quantified  by  so-called  Lyapunov  exponents.  Later  in  the  introduction 
I talk  more  about  these  exponents.  For  now  notice  that  an  A-body  system 
is  always  chaotic,  even  when  the  smooth  analogue  is  an  integrable  potential! 
However  after  a little  thought,  this  is  quite  reasonable.  An  integrable  potential 
possesses  symmetries  (as  many  as  the  degrees  of  freedom).  The  phase  space  has 
2 n dimensions  where  n is  the  degrees  of  freedom.  Because  of  the  presence  of  the 
discreteness  effects,  the  frozen- A'  analogue  does  not  possess  the  same  symmetries. 
Its  phase  space  is  more  complicated  and  it  has  nN  dimensions.  As  a result,  an 
orbit  that  is  regular  in  the  smooth  potential  can  (and  in  general  will)  be  chaotic  in 
the  frozen- A environment. 

One  might  intuitively  expect  that  as  A increases,  the  phase  space  of  a frozen- 
A system  will  resemble  more  the  phase  space  associated  with  the  smooth  potential. 
This  means  that  if  an  orbit  is  regular  in  the  smooth  potential,  then  for  very  large 
A,  in  its  frozen-A  analogue  it  should  be  regular  or  at  least  approach  regularity 
(always  in  the  sense  of  sensitive  dependence  on  initial  conditions!).  The  truth 
however  (and  the  big  surprise)  is  that  this  does  not  happen.  The  orbit  remains 
chaotic  even  for  very  large  A!  There  is  a theoretical  explanation  of  this  fact  and 
I will  talk  more  carefully  about  this  in  the  conclusions.  However  the  claim  is  that 
if  one  perturbs  an  orbit  in  an  A-body  environment,  the  time  scale  r on  which  the 
separation  of  the  perturbed  from  the  unperturbed  orbit  tends  to  grow  will  not 
diverge  (will  not  increase  without  bound)  for  A ~¥  oo.  There  is  some  disagreement 
though  as  to  whether  r should  converge  toward  an  A-independent  value  [2]  or 
whether  r should  instead  slowly  decrease  with  increasing  A [3]. 

Let’s  go  back  to  numerics.  A major  motivation  for  this  research  was  the  fact 
that  most  of  what  was  done  before  in  the  context  of  chaos  in  A-body  environments 
focused  either  on  small  A and  long  times  of  integration,  or  large  A but  short  times 
of  integration.  Little  or  no  work  had  been  done  with  large  A and  long  times  of 


integration,  a condition  needed  in  order  to  estimate  honest  Lyapunov  exponents, 
and  therefore  to  quantify  chaos  adequately. 
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A second  major  point  of  motivation  emerged  because  of  the  common  use  of  the 
a so-called  softening  parameter  (usually  denoted  as  e)  in  the  IV-body  simulations. 


softening  parameter  attractive  if  not  necessary.  If  one  looks  carefully  at  Formula 
1.3  it  is  obvious  that  if  two  particles  happen  to  come  very  close,  the  denominator 
will  become  very  large  causing  an  overflow  error  [4] . Therefore  it  is  considered  a 
proactive  behaviour  to  slightly  manipulate  the  A-body  potential  by  including  in 
the  denominator  a small  number  so  to  avoid  this  error: 


However  this  is  not  the  only  reason  the  softening  parameter  is  included.  Larger  val- 
ues of  e simply  transform  the  particle  into  a superparticle,  a group  that  is  of  many 


parameter  will  simply  put  a lower  bound  in  the  force  of  their  interaction,  which  is 
actually  what  should  happen  realistically  in  an  environment  of  A superparticles. 

The  value  of  this  idea  is  of  course  that  one  can  hope  to  simulate  a realistic  sys- 
tem of  many  bodies  using  a much  smaller  number  A in  the  numerical  simulation. 
This  is  quite  a useful  idea  but  there  is  some  consequence.  Exactly  because  the 
softening  parameter  introduces  a lower  bound  in  the  interaction  force  it  practically 
“kills”  all  the  effects  of  close  encounters  that  are  responsible  for  destabilization  of 
the  orbit  [6].  This  means  that  it  removes  all  the  chaotic  effects  that  ought  to  be 
present!  Therefore  it  is  not  a hyperbole  if  someone  views  all  simulations  done  with 
very  large  e with  a critical  eye  and  a sceptical  mind.  One  last  thing  to  mention 
about  the  softening  parameter  is  that  its  necessity  (to  avoid  overflow  errors)  makes 


First  I have  to  refer  to  the  technical  coi^iderations  that  make  the  use  of  the 


(1.4) 


particles  [5] . If  two  of  these  superparticles  happen  to  come  very  close  the  softening 


7 

it  also  necessary  to  add  one  more  element  in  the  definition  of  the  “continuum 
limit.”  The  continuum  limit  is  approached  when  e — > 0 and  N —>  oo. 

Having  described  what  is  already  known,  the  next  step  is  to  define  carefully 
what  answers  one  is  interested  in  and  what  questions  he  has  to  ask  in  order  to  find 
these  answers. 

The  main  question  is  rather  obvious:  How  close  must  a system  of  N bodies  be 
to  the  continuum  limit  (i.e.  how  large  must  N be)  in  order  for  the  dynamics  of  the 
smooth  potential  and  its  analogue  frozen- N system  to  be  indistinguishable?  The 
same  question  can  be  asked  in  different  terms:  For  what  time  scales,  i.e.  for  how 
long  the  dynamics  be  the  same?  This  is  an  interesting  question  but  very  general. 
The  term  “dynamics”  is  of  course  very  broad  and  involves  many  measures.  The 
simplest  comparison  here  would  be  in  the  context  of  the  motion  of  single  orbits, 
that  is  pointwise.  In  this  thesis  I examined  very  carefully  this  question  and  I will 
discuse  it  extensively.  In  order  to  make  it  clearer:  how  will  two  orbits,  starting 
from  the  same  initial  condition,  evolve  in  the  two  different  systems,  the  smooth 
and  the  frozen- IV?  How  do  they  diverge  in  time  and  how  does  this  divergence 
depend  on  N?  How  does  the  divergence  law  differ  for  integrable  and  nonintegrable 
potentials?  We  investigated  carefully  and  gave  answers  to  all  these  questions. 

Now  the  next  obvious  question  involves  the  quantification  of  chaos.  There  are 
two  measures  that  have  been  extensively  used  for  this  purpose.  The  first  has  been 
already  been  mentioned:  Lyapunov  exponents  [7,  8].  The  mathematical  definition 
of  this  measure  is  rather  old  and  for  completeness  and  historical  reasons  I will 
include  its  formal  definition: 


X—  lim  lim  (-)  ln( 

t-yoodz^O^t  ^5^(0)) 


(1.5) 


Of  course  numerically  this  measure  can  never  be  computed  because  of  the  limits 
(one  cannot  integrate  to  infinity).  Instead  what  was  computed  was  the  largest 
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short-time  Lyapunov  exponents,  using  a method  introduced  in  1980  by  Bennetin, 
Galgani  et.al.  [9]  which  has  since  become  classic.  This  method  simply  involves 
evolving  an  orbit  and  a second  slightly  perturbed  orbit  for  a short-time  (about  one 
dynamical  time),  then  renormalizing  the  perturbed  orbit  (to  bring  the  two  orbits 
close  again),  and  then  repeating  the  same  until  the  average  Lyapunov  exponent 
converges  to  an  almost  constant  value. 

The  second  measure  which  proved  to  be  extremely  useful  in  this  series  of 
experiments  was  the  power  spectrum  of  the  orbits  [7,  8].  It  is  well  known  that 
the  power  spectrum  of  a regular  orbit  is  characterized  by  a few  almost  discrete 
frequencies  (theoretically  they  are  discrete  but  this  can  never  be  the  case  for 
numerical  experiments,  unless  one  integrates  for  t -+  oo,  a numerical  impossibility 
of  course).  On  the  other  hand  chaotic  orbits  are  characterized  by  a much  more 
complicated  spectrum,  that  may  include  a few  ‘main’  frequencies,  but  also  includes 
lots  of  other  frequencies  with  smaller  but  collectively  significant  power.  These  extra 
frequencies  give  some  “grassy”  look  to  the  spectrum. 

There  are  still  two  questions  of  importance.  The  first  involves  the  concept 
of  “chaotic  mixing”  (i.e.  the  exponential  divergence  of  initially  localised  chaotic 
orbits)  [10,  11,  12,  13,  14].  How  does  chaotic  mixing  take  place  in  an  IV-body 
environment?  How  does  it  depend  on  N?  And  what  are  the  answers  to  the  same 
questions  for  regular  orbits?  An  extensive  investigation  of  these  issues  was  also 
undertaken. 

The  last  question  is  one  of  major  practical  and  potentially  theoretical  interest. 
Is  it  possible  to  simulate  motion  in  a frozen- N environment  using  a Langevin 
equation  [15,  16]?  This  means  that  one  has  to  include  in  the  canonical  equations 
of  motion  some  term  that  represents  a Gaussian  process  with  a white  noise  power 
spectrum  in  order  to  simulate  the  random  kicks  a test  particle  “feels”  as  it  moves 
in  the  frozen- N environment.  Could  such  an  approach  give  results  similar  to  the 
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results  of  real  iV-body  simulations?  A number  of  experiments  were  also  performed 
here  with  extremely  interesting  results. 


CHAPTER  2 

MATHEMATICAL  FORMULAS 


Hitherto  I have  refered  to  what  information  already  exists  in  work  that  has 
been  done  by  other  researchers.  I have  also  defined  carefully  what  the  questions  are 
and  I have  given  a broad  outline  about  how  to  get  these  answers,  what  experimen- 
tal procedures,  measures,  and  criteria  are  necessary.  Next  I will  get  more  specific 
about  the  experimental  processes  used  for  this  thesis. 

Firstly  it  is  necessary  to  give  in  detail  all  the  mathematical  formulas  that 
describe  the  smooth  potentials  which  were  chosen  to  be  used  as  well  as  brief 
information  pertaining  to  the  reasons  those  potentials  were  chosen. 

In  this  investigation  four  smooth  potentials  were  used  for  the  experiments. 
They  are  the  following: 


2.1  Plummer  Potential 

This  is  a potential  as  old  as  well-known.  It  was  first  used  in  1911  by  H.  C. 
Plummer  to  fit  observations  of  globular  clusters  [5].  This  potential  is  integrable. 
Therefore  all  orbits  are  regular.  Mathematically  it  is  expressed  as: 


$P(r)  = 


GM 

-\/r2  + b2 


(2.1) 


This  potential  can  be  generated  via  Poisson’s  equation  by  the  following  density 
function: 


Pp(r) 


(2.2) 
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with  G the  gravitational  constant,  M the  mass  of  the  system  and  b a core  radius 
parameter.  For  the  experiments  performed  units  were  chosen  so  that:  G = M = 
b=  1. 


2.2  Three-dimensional  harmonic  oscillator 
This  is  also  an  integrable  potential,  and  it  is  extremely  simple  and  computa- 
tionally fast.  Its  mathematical  formula  is: 

<FE(r)  = $0  + \ Kx2  + uly2  + u2cz2)  (2.3) 


where  ua,  ujb,  and  uc  are  the  frequencies  in  the  three  different  directions  and  4>0  is 
constant.  The  values  chosen  for  the  experiments  performed  were:  4>0  ~ — 1-00608, 
id a ~ 0.4663,  ojb  ~ 0.5508,  and  ujc  « 0.6753. 

The  mass  density  associated  with  this  potential  is: 


Pe( r) 


3 M 
Airabc 


x m2 


(2.4) 


with 


m2  = 


x2  y 2 z2' 

a2  b2  c2  , 


(2.5) 


with  frequencies  ua,  Ub,  ujc,  related  to  the  axes  values  a,  b , c,  in  terms  of  incomplete 
elliptic  integrals  [30].  The  values  we  chose  for  of  a,  6,  and  c,  were  a = 1.95, 
b = 1.50,  and  c = 1.05. 


2.3  Three-dimensional  harmonic  oscillator  plus  a perturbation  factor 
The  perturbation  factor  simulates  a black  hole  at  the  center  of  the  system. 
Because  of  the  perturbation  this  potential  is  nonintegrable  and  admits  large  regions 
of  global  stochasticity.  It  has  been  found  that  the  vast  majority  of  orbits  (restricted 
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energetically  to  m < 1)  in  this  potential  are  chaotic  [31]: 

GMbh 


Qrh  = — 


Vr2  + 


=2’ 


(2.6) 


with  Mbh  the  mass  of  the  black  hole  and  e a softening  parameter  introduced  to 
avoid  overflow  errors.  G represents  as  usual  the  gravitational  constant  and  for 
the  choice  of  the  units  used  its  value  is  G = 1.  The  choices  for  the  perturbation 
parameters  were:  e = 10“3,  and  Mbh  = 10 _15M  ~ 0.0316228M.  The  parameters 
associated  with  the  integrable  part  were  the  same  as  in  the  previous  potential. 

The  mass  density  associated  with  this  potential  is: 


Pe{  r)  = Pe  + 


GM, 


BH 


—3  r2 


47T 


+ 


1 


(r2  + e2)a  (r2  + e2)2 


(2.7) 


2.4  A triaxial  generalisation  of  Dehnen  potential 
This  potential  does  not  have  an  analytic  form,  but  its  values  can  be  computed 
numerically  via  Poisson’s  equation.  The  density  associated  with  this  potential  is 
[32]: 


PD  (r) 


47t  abc 


(1  + m) 


(4-7) 


(2.8) 


with 


m2  = 


x2  y 2 

^ + — + 


b2 


(2.9) 


where  7 and  a,  b,  and  c are  free  parameters.  The  values  chosen  for  the  experiments 
were:  7 = 1,  a = 1.0,  b = ^5. 0/8.0,  and  c = 0.5.  These  values  have  been  used  to 
simulate  a triaxial  galaxy  with  a central  cusp  [33] . It  has  been  found  numerically 
that  for  these  values  this  potential  admits  a coexistence  of  large  measure  of  both 
regular  and  chaotic  orbits[33,  34]. 

For  completeness  I will  give  again  here  the  mathematical  formulas  associated 
with  the  potential  and  the  mass  density  of  frozen- iV  systems  generated  by  sampling 
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the  mass  densities  of  smooth  potentials.  The  potential  is: 

t / \ GM v ^ i 

<Mr)  = — 2^ 


(2.10) 


N \/(r  — ri)2  + e2 

(where  I have  included  the  softenting  parameter)  and  it  is  generated  by  the  density 


of  N discrete  masses: 


p(r)  = ~r‘)  (2-11) 

t=i 

The  main  experiment  performed  was  quite  straightforward:  Integrate  some  initial 
condition(s)  (test  particles)  in  a smooth  potential.  Generate  frozen- ./V  distributions 
that  correspond  to  this  potential.  Then  integrate  the  same  initial  condition(s)  in 
the  frozen- TV  environment.  For  every  intergation  (smooth  or  frozen- TV)  two  kinds  of 
data  were  recorded:  (a)  the  orbital  data,  and  (b)  the  largest  short-time  Lyapunov 
exponents  associated  with  the  orbits  of  the  test  particles. 

The  integrations  for  the  smooth  potentials  were  performed  by  an  adaptive 
time  step  Runge-Kutta  algorithm  [40] . In  the  case  of  the  Dehnen  potential  a code 
written  by  C.  Siopis  for  his  Ph.D  dissertation  was  used  [35].  The  integrations  in  the 
frozen- N environments  were  realised  by  a direct  summation  code  written  by  the 
author.  All  integrations  conserved  energy  to  at  least  one  part  in  10-4.  The  routine 
that  computed  the  Lyapunov  exponents  was  based  on  a well  known  algorithm  by 
Benettin,  Galgani  et.al  [9].  The  sampling  of  the  density  distributions  was  realized 
for  the  Plummer  potential  by  a von  Neumann  rejection  method,  and  for  the  rest  of 
the  cases  by  a code  written  by  the  author  that  was  based  on  a new  algorithm. 


CHAPTER  3 

INTEGRABLE  POTENTIALS 


Since  the  mathematical  formulation  of  the  numerical  experiments  has  been  set 
up  we  are  ready  to  proceed  to  the  description  of  the  experiments.  In  this  chapter 
I will  talk  exclusively  about  the  experiments  and  results  related  to  integrable 
potentials.  In  the  next  chapter  I will  proceed  to  the  nonintegrable  potentials. 

3.1  Divergence  Laws 

The  first  series  of  experiments  were  performed  using  the  Plummer  potential. 
The  purpose  was  to  identify  first  the  dynamical  behaviour  in  a generic  integrable 
system,  without  the  complications  introduced  by  the  chaotic  behaviour  of  noninte- 
grable potentials.  The  second  step  was  to  perform  the  same  experiments  in  the  also 
integrable  (but  nongeneric)  three-dimensional  harmonic  oscillator  and  see  if  similar 
behaviour  was  observed.  Then  having  a clear  picture  about  what  happens  in  the 
integrable  case  we  could  proceed  to  the  more  complicated  nonintegrable  case. 

All  integrations  were  performed  for  evolution  times  that  physically  correspond 
to  about  100  dynamical  times.  Such  long  times  were  both  necessary  and  important 
for  the  fair  treatment  of  the  problem  because:  (a)  100 tD  (for  galaxies)  is  a time 
scale  equivalent  to  the  age  of  the  universe,  so  it  is  the  natural  time  scale  for 
the  evolution  we  wanted  to  simulate,  and  (b)  (a  more  technical  consideration) 
integration  times  of  such  length  are  required  numerically  since  they  are  necessary 
for  the  Lyapunov  exponents  to  converge. 

The  obvious  first  comparison  was  between  the  evolution  of  initial  conditions 
(r0,  v0)  integrated  in  the  smooth  potential  and  in  the  corresponding  frozen- IV 
environment.  This  comparison  should  be  performed  for  different  N and  for  variable 
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values  of  the  softening  parameter.  This  way  a clear  first  picture  would  emerge 
about  the  validity  of  the  continuum  limit.  Twenty  different  frozen- N distributions 
were  generated  for  each  choice  of  N.  The  number  of  particles  for  the  frozen- iV 
environments  varied  from  N = 102  5 to  N = 105  5.  The  softening  parameter  also 
varied  from  e = 10_1  to  e = 10-5.  The  closer  to  zero  e is,  the  more  “honest”  the 
frozen- TV  simulation  is  considered. 

Intuitively  one  would  expect  that  for  sufficiently  short  times  the  two  orbits 
will  coincide,  but  that  for  longer  times  they  should  diverge  because  the  smooth  and 
the  frozen- iV  potential  are  only  approximately  similar.  The  differences  the  body 
experiences  during  his  motion  in  the  two  environments  should  eventually  translate 
into  a divergence  of  the  two  orbits. 

Therefore  two  meaningful  questions  can  easily  arise:  (a)  How  do  the  position 
and  velocity  divergence  between  the  smooth  characteristic  and  the  frozen- TV  orbit 
evolve  in  time?  Do  they  evolve  as  a power  law  or  exponentially?  (b)  How  does 
the  growth  time  of  this  divergence  depend  on  the  number  of  particles  N ? The 
answer  to  these  questions  is  quite  critical  for  the  validity  of  the  continuum  limit 
at  a pointwise  level.  Equally  important,  however,  is  to  ask  the  same  two  questions 
about  the  angular  momentum,  which  is  an  invariant  of  motion  for  the  Plummer 
potential  since  it  is  a spherically  symmetric  system. 

First  let’s  consider  these  questions  at  the  most  elementary  level  and  examine 
a crude  comparison  between  two  orbits,  starting  from  the  same  initial  condition, 
but  with  the  first  evolved  in  the  smooth  potential  and  the  second  in  the  frozen- iV 
environment. 

It  is  probable  that  the  first  naive  expectation  (or  guess!)  someone  would  have 
about  the  law  that  governs  the  divergence  of  the  two  orbits  in  the  two  different 
regimes  would  be  that  it  is  exponential,  since  it  is  known  that  the  frozen- A?”  orbits 
are  chaotic.  However  exponential  divergence  was  definitely  not  observed.  Rather 
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t/to  t/t0 

Figure  3.1:  (a)  The  deviation  A r — |r  — rs|  from  the  smooth  potential  characteristic 
for  one  representative  frozen-A  orbit  evolved  with  e = 10-4  and  A = 103.  (b)  An- 
other frozen- A'  orbit  evolved  with  the  same  e and  A for  the  same  initial  condition, 
(c)  and  (d)  the  same  for  A = 105. 
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what  was  observed  was  a power  law  in  time,  a striking  result!  This  was  true  not 
only  for  small  values  of  e,  where  frozen- N orbits  are  very  chaotic,  but  also  for  larger 
values  of  e,  where  orbits  tend  to  be  more  regular  because  of  the  suppression  of 
chaos.  Figure  3.1  exhibits  the  quantity  A r — |r  — rs|  for  a typical  orbit  evolved  in 
the  frozen-iV  environment  with  e = 10-4.  The  first  two  panels  are  for  N = 103,  and 
the  latter  two  for  N = 105.  It  is  evident  that,  rather  than  exponential,  the  growth 
of  A r is  roughly  linear  in  time. 

This  divergence  continues  until  the  typical  separation  has  become  comparable 
to  the  size  of  the  configuration  space  region  to  which  the  orbits  are  confined,  and 
the  frozen-iV  orbit  has  become  completely  decorrelated  in  appearance  from  the 
smooth  potential  characteristic.  This  means  that  the  frozen-iV  orbits  do  not  stay 
close  to  the  smooth  characteristic,  which  would  be  the  case  if  the  divergence  was 
terminated  at  comparatively  small  separations. 

Now  let’s  examine  the  full  statistics  of  the  problem  using  moments  that  involve 
all  20  simulations  in  the  frozen- N configurations.  Two  types  of  moments  were  used. 
First  let’s  denote  the  positions  and  velocities  of  the  n — 20  different  trajectories 
(for  one  initial  condition)  in  the  20  different  frozen-iV  environments  as  (r i(t),v;(£)), 
{i  — 1,  ...,n),  and  the  position  and  velocity  of  the  smooth  characteristic  (for  the 
same  initial  condition)  as  (r s(t),v3(t)).  The  first  type  of  moment  that  was  used  (for 
the  position  variables)  included  the  following: 

= (3-1) 

i=l 

and 

Dr 2 = <|r*  - (r)|2)  = ^ ^ |r,  - (r)|2 


(3.2) 
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The  second  type  included: 

Ar2  = (|r-r,|2)  = -V  |r*  - r,|2  (3.3) 

n z— ■ ' 
i= 1 

and 

<5r2  = (|(r)  - rs|2)  (3.4) 

The  same  quantities  were  also  generated  for  the  velocity  v and  for  the  angular 
momentum  J = rxv.  It  was  found  that  the  quantities  Dr,  Ar,  and  Sr  exhibit 
comparatively  similar  evolutions.  For  this  reason  only  Ar  will  be  discussed. 

Before  I proceed,  I will  refer  very  briefly  to  an  analytical  consideration.  Let’s 
denote  by  D2  the  mean  value  of  r2  associated  with  the  smooth  characteristic.  Now 
given  that  the  frozen- TV  orbits  conserve  energy,  the  average  distance  from  the  origin 
should  be  the  same  as  for  the  smooth  characteristic,  (that  is  D2)  so  that: 

(r2(f))  — > R2s  for  t — > oo  (3.5) 

Also  one  expects  that: 

(r (t)  ■ r s(t))  -4  0 for  t — » oo  (3.6) 

so  that  it  may  be  infered  that  Ar2  -»  2D2. 

It  follows  that  one  can  use  2D2  as  a normalization  factor  and  plot  the  quantity 
Ar/ y/2 R2S  instead  of  Ar.  The  normalized  quantity  then  should  saturate  towards 
unity,  reflecting  the  fact  that  the  orbits  have  become  completely  different  from  each 
other. 

This  is  what  can  actually  be  seen  in  Figure  3.2.  The  left-hand  two  panels 
show  the  evolution  for  integrations  with  e = 10-4  for  N = 103  and  N = 105.  The 
right-hand  panels  exhibit  the  same  data  but  now  subjected  to  a boxcar  averaging. 


Ar/(2'/!Rs)  Ar/(2'/!Rs) 
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Figure  3.2:  (a)  The  quantity  Ar/y/2K*  for  frozen- N simulations  with  N — 103  and 
e = 10-4.  (b)  The  same  data  subjected  to  boxcar  averaging  over  an  interval  t = tp. 
(c)  and  (d)  The  same  for  N = 105. 
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Figure  3.3:  Best  fit  values  of  the  time  scale  tG(N)  associated  with  the  di- 
vergence of  smooth  and  frozen- iV  orbits  for  two  different  initial  conditions: 

A r/Ra  = Av/Vs  = t/tG-  The  dashed  line  has  slope  1/2,  corresponding  to  an  TV1/2 
dependence. 

Let’s  discuss  about  what  happens  on  these  plots.  I have  already  explained  why 
the  saturation  should  happen  at  a value  of  unity.  The  large  envelopes  associated 
with  the  curves  in  the  left-hand  panels  reflect  the  fact  that  at  late  times  the  frozen- 
N orbits  oscillate  about  the  value  of  unity.  Now  two  additional  facts  are  quite 
obvious  from  this  graph.  The  first  is  that  the  evolution  in  time  is  linear.  And  the 
second  is  that  the  growth  time  of  the  divergence  (i.e.  the  slope  of  the  curve  from 
time  zero  to  the  saturation  time)  depends  on  the  number  N of  the  frozen  particles. 
It  is  obvious  that  the  divergence  grows  more  slowly  as  N increases,  which  is  what 
someone  would  intuitively  expect,  since  larger  N means  that  the  frozen- A potential 
resembles  the  smooth  potential  more  closely.  Figure  3.3  exhibits  \og10(tG/tD)  as  a 
function  of  log10  N for  two  different  initial  conditions  integrated  for  e = 10~5. 
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Thus  the  mathematical  formula  that  expresses  the  observed  behaviour  can  be 
written  as: 


A r 
r. 


t 

tG 


(3.7) 


Here  to  represents  the  growth  time  and  is  given  by: 

to  ~ AcNptD 

where  A is  of  order  unity  and  p « 0.5.  The  complete  law  is  thus: 

Ar  1 t 
— oc 


(3.8) 


The  same  law  applies  for  Av: 


y/NtD 


Av  1 t 
— oc 
v„ 


(3.9) 


(3.10) 


VNtD 

The  square  root  dependence  on  N of  both  Ar  and  Av  might  suggest  diffusion  but 
the  time  dependence  is  definitely  not  square  root.  However  such  a dependence 
in  time  was  observed  for  the  corresponding  statistical  quantities  associated  with 
angular  momentum.  There  (for  small  e)  results  suggest  the  following  dependence: 


A J2 


t 


J2s  tj 

where  Js  is  the  magnitude  of  the  angular  momentum  associated  with  the  smooth 
characteristic  and 


(3.11) 


tj  ~ AjNptc>  (3.12) 

where  Aj  ~ 3 AG  and  p « 0.5.  The  dependence  of  \ogw(tj/tD)  versus  log10  N can 
be  seen  in  Figure  3.4. 
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Figure  3.4:  Best  fit  values  of  the  time  scale  tj(N)  associated  with  changes 
in  angular  momentum  for  frozen- iV  orbits  for  two  different  initial  conditions: 
AJ2/J2  = t/tj.  The  dashed  line  has  slope  1/2,  corresponding  to  an  N1/2  depen- 
dence. 
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3.2  Average  Short-time  Lyapunov  Exponents 

Largest  short-time  Lyapunov  exponents  have  been  used  extensively  as  the 
principal  diagnostic  of  chaos.  The  aim  here  was  to  compute  the  average  short-time 
Lyapunov  exponent  characterizing  orbits  starting  from  the  same  initial  condition 
and  evolved  in  the  20  different  frozen- N configurations,  for  different  values  of  N 
and  e.  The  dependence  of  the  average  short-time  Lyapunov  exponent  on  N and  e 
can  be  seen  in  Figures  3.5  and  3.6.  It  is  obvious  from  Figure  3.5  that  the  size  of  the 
Lyapunov  exponent  decreases  as  e increases.  This  is  an  expected  behaviour.  Since 
the  potential  is  integrable  the  observed  chaos  is  generated  completely  by  the  close 
encounters  of  the  test  particles  with  the  field  particles.  By  increasing  e,  a minimum 
impact  parameter  is  imposed  to  the  system,  so  that  the  destabilization  of  the  orbit 
(associated  with  the  close  encounters)  that  leads  to  chaotic  behaviour,  is  suppresed. 

However  when  e is  sufficiently  small,  its  precise  value  seems  to  be  immaterial. 
This  happens  because  as  long  as  e is  smaller  than  the  value  of  the  closest  separa- 
tion between  the  test  particle  and  the  field  particles  during  the  integration,  the  test 
particle  feels  practically  an  unsoftened  potential.  The  obvious  point  derived  from 
these  results  is  that  the  introduction  of  large  softening  definitely  suppresses  the 
perturbations  that  can  lead  to  chaotic  behaviour. 

When  e is  large  the  value  of  the  Lyapunov  exponent  decreases  with  increasing 
N . The  reason  is  obvious:  For  a system  of  fixed  size,  as  N increases  the  typical 
distance  between  particles  decreases.  For  very  small  N this  distance  is  still  large, 
larger  than  the  value  of  the  softening  parameter,  so  that  the  close  encounters 
that  generate  chaos  are  not  suppressed.  However  as  N increases,  the  interparticle 
distance  decreases.  When  this  distance  becomes  smaller  than  e,  many  weak 
encounters  get  suppresed,  a fact  that  causes  the  observed  decrease  in  the  value  of 
the  Lyapunov  exponents. 
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Figure  3.5:  Mean  short-time  Lyapunov  exponent  (x)  as  a function  of  softening  pa- 
rameter e for  N = 105  (solid  line),  N — 1045  (dotted),  N = 104  (dashed),  N = 103-5 
(dot-dashed),  N = 1030  (triple-dot  dashed),  and  N = 102  5 (broad  dashed).  The 
integrations  were  all  performed  for  a single  ‘typical’  initial  condition. 
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Figure  3.6:  Mean  short-time  Lyapunov  exponent  (x)  as  a function  of  particle  num- 
ber N for  e = 10-5  (solid  line),  e = 10-4  (dotted),  e = 10-3  (thin-dashed),  e = 10-2 
(dot-dashed),  and  e = 10_1  (triple-dot  dashed),  all  computed  for  the  initial  condi- 
tion used  to  generate  Figure  3.5.  The  short-time  (x)  for  a different  initial  condition 
corresponding  to  a smooth  radial  orbit,  again  evolved  with  e = 10~5,  is  indicated  by 
the  curve  with  thick  dashes. 
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However  when  e is  very  small  the  Lyapunov  exponents  seem  almost  indepen- 
dent of  N.  It  is  possible  to  understand  qualitatively  the  physical  reason  for  this. 
Integrable  potentials  do  not  exhibit  chaos.  This  means  that  all  the  chaos  observed 
in  their  frozen- iV  analogues  must  be  associated  with  discreteness  effects  that 
perturb  the  orbits  because  of  close  encounters.  This  suggests  that  the  time  scale 
for  the  growth  of  a small  initial  perturbation  can  be  calculated  if  we  consider  the 
tidal  effects  associated  with  two  particles  separated  by  some  typical  interparticle 
distance.  The  tidal  acceleration  is: 

(...  , r Gw 

Sr  = (Sr  ■ V)a  ~ — — 5r  (3.13) 

where  r is  the  separation  and  m the  particle  mass.  Suppose  however  that: 

r ~ n~3  ~ N~^ Rsys  (3.14) 


with  n the  characteristic  number  density  and  Rsys  the  size  of  the  system.  It  then 
follows  dimensionally  that  the  time  scale  that  characterizes  the  interaction  is: 


U 


r>sj 


1 

7&P 


(3.15) 


which  is  independent  of  N and  practically  comparable  to  the  dynamical  time  (also 
independent  of  N). 


3.3  Power  Spectra  and  Complexities 

The  results  described  in  the  previous  paragraph  about  the  average  short-time 
Lyapunov  exponents  suggest  clearly  that  for  small  e the  orbits  remain  chaotic 
regardless  of  the  value  of  N.  Let’s  now  have  a look  at  the  visual  appearance  of  a 
typical  orbit  in  the  frozen-iV  configuration  for  different  values  of  N (102  5 — 105'5) 
and  e = 10“5  (Figure  3.7).  The  final  panel  exhibits  the  behaviour  of  the  orbit 
in  the  smooth  potential  environment.  It  is  very  clear  from  this  picture  that  as  N 
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increases:  (a)  the  orbit  looks  less  chaotic,  and  (b)  it  resembles  more  closely  the 
smooth  characteristic. 

In  order  to  visualize  in  a more  quantitative  fashion  this  transformation 
towards  smoother  orbital  behaviour  as  N increases,  one  can  consider  the  power 
spectrum  of  the  orbit.  The  power  spectra  that  correspond  to  the  orbits  in  Figure 
3.7  can  be  seen  in  Figure  3.8.  First  let’s  have  a look  at  the  power  spectrum  of  the 
smooth  orbit.  It  is  clear  that  the  power  is  concentrated  into  two  main  frequencies 
that  are  very  sharp  (this  is  expected  for  a regular  orbit).  Now  for  very  small  N 
the  power  spectrum  exhibits  a variety  of  different  frequencies,  something  that 
characterizes  a chaotic  orbit.  However,  as  N increases  the  power  spectrum  changes 
to  resemble  more  closely  the  spectrum  of  the  smooth  characteristic,  with  the  power 
concentrated  into  smaller  number  of  frequencies,  and  the  “grass”  of  frequencies 
associated  with  chaos  largely  eliminated. 

In  order  to  quantify  this  result  more  objectively,  the  next  logical  step  is  to 
compute  the  complexity  of  the  orbits,  i.e.  the  number  of  frequencies  in  the  Fourier 
spectrum  that  contains  an  appreciable  amount  of  power  [29] . Two  measures  of 
complexity  were  defined.  The  first  /0.i  is  the  sum  of  the  number  of  frequencies 
/o.i,x)  /o.i.j/,  and  /o.i,z  which  have  more  than  10%  as  much  power  as  the  peak 
frequencies  for  wx,  cuy,  and  a iz\ 


The  second  measure  k0, 95  is  defined  correspondingly  for  the  sum  of  the  number  of 
frequencies  required  to  capture  95%  of  the  power  in  the  x,y,  and  z directions: 


/o.l  — fo.l,x  + Zb. l,t/  + fo.i,z 


(3.16) 


(3.17) 


The  graph  of  both  measures  versus  logwN  can  be  seen  in  Figure  3.9.  With 
increasing  N both  these  quantities  decrease. 
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Figure  3.7:  The  x-y  projection  of  representative  frozen- iV  orbits  generated  from 
the  same  initial  condition,  evolved  for  t = 25 tD  with  e = 10~5.  (a)  N = 102  5. 
(b)  N = 103.  (c)  N = 103-5.  (d)  N = 104.  (e)  N = 1045.  (f)  N = 105.  (g) 
N = 105-5.  (h)  The  x-y  projection  of  the  same  initial  condition  evolved  in  the 

smooth  potential. 
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Figure  3.8:  (a)  The  Fourier  transformed  |x(cu)|  for  one  frozen- N integration  of  the 
initial  condition  used  to  generate  Figure  3.7,  evolved  with  e = 10-5  and  N = 102'5. 
(b)  The  same  for  N = 103.  (c)  N = 103'5.  (d)  N — 104.  (e)  N — 104  5.  (f)  N = 105. 
(g)  N = 105'5.  (h)  |a:(u;)|  for  a characteristic  in  the  smooth  potential  with  the 

same  initial  condition,  with  data  recorded  at  the  same  intervals  for  the  same  total 
integration  time. 
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Figure  3.9:  Two  probes  of  the  complexity  of  frozen- N orbits  for  an  ensemble  of  or- 
bits with  the  same  initial  condition  evolved  with  e = 10-5.  The  solid  curve  exhibits 
/o.i,  the  number  of  frequencies  which  have  power  equal  to  at  least  10%  of  the  power 
in  the  peak  frequencies.  The  dashed  curve  exhibits  /co.95,  the  number  of  frequencies 
required  to  capture  95%  of  the  total  power.  The  vertical  lines  show  /0. 1 and  fco.95 
for  a smooth  characteristic  generated  identically  from  the  same  initial  condition, 
thus  exhibiting  the  intrinsic  limitations  associated  with  the  discrete  time  series  of 
data  points. 


CHAPTER  4 

NONINTEGRABLE  POTENTIALS 
4.1  Divergence  Laws 

In  order  to  investigate  the  law  characterizing  the  divergence  of  a chaotic  orbit 
moving  in  a smooth  nonintegrable  potential  from  an  orbit  starting  from  the  same 
initial  condition  but  evolving  in  the  frozen-N  analogue  of  the  smooth  potential,  the 
three-dimensional  harmonic  oscillator  plus  perturbation  potential  (Equation  2.6) 
was  used. 

Here  the  same  statistical  moments  were  used  as  for  the  integrable  potential 
case.  The  definitions  can  be  seen  in  the  Chapter  3 (Equations  3. 1-3.4)  and  there 
is  no  need  to  repeat  them  here.  A number  of  simulations  were  performed  for  a 
statistical  sample  of  800  test  particles  starting  from  initial  conditions  that  in 
the  smooth  potentials  evolved  as  “wildly”  chaotic  orbits  (the  typical  Lyapunov 
exponent  was  x — 0.055).  The  results  of  these  simulations  are  exhibited  in  Figure 
4.1.  It  is  obvious  that  the  position  divergence  evolves  linearly  in  time,  the  same 
way  that  it  evolved  for  the  integrable  case.  The  same  is  true  for  the  velocity 
divergence. 

However  there  is  an  important  difference  between  the  integrable  and  noninte- 
grable cases.  The  growth  time  of  the  divergence  was  found  in  the  integrable  case  to 
depend  on  the  square  root  of  N.  In  the  nonintegrable  case  this  is  no  longer  true. 
Now  the  dependence  is  the  following: 

tG  « AG(\nN)tD  (4.1) 
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Figure  4.1:  The  mean  separation,  Dr , between  frozen-iV  orbits  and  smooth  char- 
acteristics with  the  same  initial  conditions,  computed  for  ensembles  of  800  chaotic 
initial  conditions  evolved  in  the  potential  (2.6)  for  variable  N.  (a)  N = 102'5.  (b) 
N = 103.  (c)  N = 103'5.  (d)  N = 104.  (e)  N = 1045.  (f)  N = 105.  (g)  N = 1055. 
(h)  The  growth  rate  ta{N)  extracted  from  the  preceding  panels,  assuming  a linear 
growth  law.  The  solid  line  overlays  a least  squares  fit  = A + B log10  N. 


where  Ag  is  of  order  of  unity.  So  the  complete  divergence  law  for  the  case  of 
nonintegrable  potentials  can  be  expressed  as: 
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A r 1 t 

rs  ^ In  N tD 

and  for  the  velocity: 

Av  1 t 

— oc 

vs  In  N tp 


(4.2) 


(4.3) 


where  (r s(t),vs(t))  express  the  position  and  velocity  of  the  smooth  characteristic. 

The  last  panel  of  Figure  4.1  shows  the  logarithmic  dependence  of  the  growth 
time  tg  on  N. 


4.2  Average  Short-time  Lyapunov  Exponents 
Average  short-time  Lyapunov  exponents  were  computed  for  orbits  evolved  in 
the  frozen- A environments  associated  with  the  nonintegrable  three-dimensional 
harmonic  oscillator  plus  a perturbation  potential.  The  black  hole  mass  was 
MBh  = 10“ 15 , i.e.  substantially  smaller  than  unity.  Initial  conditions  for  both 
wildly  chaotic  (orbits  that  explore  all  of  the  phase  space  energetically  available 
to  them)  and  “sticky”  chaotic  orbits  (orbits  that  remain  for  some  time  in  a small 
region  of  the  phase  space  but  eventually  escape  to  explore  all  the  available  phase 
space  energetically  available  to  them)  were  integrated.  Also  initial  conditions  that 
correspond  to  regular  orbits  in  the  smooth  potential  regime  were  integrated  in 
frozen- iV  environments  associated  with  the  three-dimensional  harmonic  oscillator 
without  the  black  hole,  a completely  integrable  potential.  The  average  Lyapunov 
exponents  for  four  different  initial  conditions  for  every  choice  of  N can  be  seen  in 
Figure  4.2.  Several  very  important  points  can  be  derived  from  this  picture,  (a) 
The  sizes  of  the  Lyapunov  exponents  for  the  integrable  and  nonintegrable  case 
are  comparable.  This  means  that  someone  knowing  only  the  size  of  the  Lyapunov 
exponent  of  an  orbit  that  has  been  evolved  in  a frozen- N environment  cannot 
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Figure  4.2:  Estimates  of  the  largest  Lyapunov  exponent  for  orbits  evolved 
in  frozen- N realisations  of  a homogeneous  ellipsoid:  integrable  initial  condi- 
tions evolved  with  Mbh  = 0 (solid  line),  ‘sticky’  initial  conditions  evolved 

with  Mbh  = 10-1,5  (dashed  line),  and  ‘wildly  chaotic’  initial  conditions  with 

Mbh  — 10-1'5  (dot-dashed  line). 

know  if  this  orbit  will  evolve  regularly  or  chaotically  in  the  smooth  potential 
environment,  (b)  For  different  numbers  of  particles  N the  size  of  the  Lyapunov 
exponents  for  nonintegrable  potentials  is  again  nearby  constant.  Recall  that 
the  same  is  true  for  integrable  potentials  (both  the  three-dimensional  harmonic 
oscillator  and  the  Plummer  potential  in  the  last  chapter  exhibited  the  same 
behaviour),  (c)  The  typical  values  of  the  Lyapunov  exponents  of  orbits  evolved 
in  the  frozen- N configurations  (~  0.8)  are  very  large  compared  to  the  values 
associated  with  typical  chaotic  orbits  evolved  in  the  smooth  potential  (~  0.055  for 
wildly  chaotic  orbits  and  ~ 0.022  for  “sticky”  chaotic  orbits). 

The  last  point  is  very  important.  The  Lyapunov  exponents  that  characterize 
the  chaotic  orbits  in  smooth  environments  are  associated  with  the  effect  the  bulk 
potential  has  on  the  orbit.  Therefore  they  probe  a macroscopic  effect.  On  the 
other  hand  the  large  Lyapunov  exponents  that  characterize  the  orbits  in  frozen- TV 
environments  are  associated  with  the  discreteness  effects,  that  is,  a microscopic 
effect.  The  chaotic  effect  associated  with  the  bulk  potential  is  still  present  in  the 
frozen- N configuration  but  it  is  so  small  that  it  is  practically  “washed  out”  by  the 
microscopic  chaos  generated  by  the  discreteness  effects. 


35 


4.3  Power  Spectra  and  Complexities 
The  same  effects  that  were  observed  for  the  power  spectra  and  visual  ap- 
pearance for  the  integrable  potentials,  were  also  observed  for  the  nonintegrable 
potentials.  A chaotic  orbit  evolved  in  a smooth  potential  is  characterized  by  a 
very  complex  power  spectrum.  Therefore  the  motion  is  characterized  by  many 
frequencies  (the  complexity  is  large).  The  same  orbit  evolved  in  the  frozen- N ana- 
logue of  the  nonintegrable  smooth  potential  will  have  an  even  more  complex  power 
spectrum  if  N is  small  but  as  N increases  and  the  continuum  limit  is  approached 
this  power  spectrum  will  resemble  closer  the  spectrum  of  the  orbit  evolved  in  the 
smooth  potential.  The  same  of  course  happens  to  the  visual  appearance  of  the 
orbit. 


4.4  Transitions  Between  Regular  and  Chaotic  Behaviour 
In  a frozen- TV  environment  it  is  possible  for  an  orbit  to  “jump”  from  the 
regular  regime  to  the  chaotic  one  and  vise  versa.  This  happens  because  discretenss 
effects  can  push  the  orbit  from  a path  that  corresponds  in  the  smooth  regime  to 
a chaotic  orbit,  to  a new  path  that  corresponds  in  the  smooth  regime  to  a regular 
orbit,  or  vice  versa.  This  of  course  is  a behaviour  that  can  be  observed  only  for 
cases  of  potentials  that  admit  significant  measures  of  both  regular  and  chaotic 
orbits.  One  example  of  such  a potential  is  the  triaxial  generalisation  of  the  Dehnen 
potential  described  in  Chapter  2. 

The  first  hint  that  such  transitions  do  happen  is  given  by  the  time  evolution  of 
the  divergence  Dr  for  two  different  ensembles  of  initial  conditions,  one  that  evolves 
regularly  and  one  that  evolves  chaotically  in  the  smooth  regime.  The  evolution  of 
the  divergence  can  be  seen  in  Figure  4.3.  What  is  clear  is  that  the  growth  time  for 
both  ensembles  is  comparable  especially  for  small  N.  For  larger  N the  divergence 
grows  slower  for  regular  orbits  as  it  is  expected.  The  fact  that  for  smaller  N the 
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Figure  4.3:  The  quantity  Dr  for  two  different  initially  localised  ensembles  of 
‘sticky’  chaotic  orbits  with  the  same  (low)  energies,  each  evolved  in  frozen-iV  re- 
alisations of  the  Dehnen  potential,  (a)  A ‘regular’  ensemble  with  N — 104.  (b)  A 

‘chaotic’  ensemble  with  N = 104.  (c)  The  ‘regular’  ensemble  with  N = 104,5.  (d) 
The  ‘chaotic’  ensemble  with  N = 104,5.  (a)  The  ‘regular’  ensemble  with  N = 105. 
(b)  The  ‘chaotic’  ensemble  with  N — 105. 


growth  times  are  comparable  can  be  explained  as  transitions  of  significant  numbers 
of  orbits  from  regular  to  chaotic  and  vice  versa. 

In  order  to  quantify  more  carefully  these  transitions  the  following  experiment 
was  performed.  Snapshots  of  the  orbital  data  of  the  orbits  evolved  in  frozen- N 
environments  (i.e.  the  positions  and  velocities  of  the  particles  at  times  to,  t\,  t2 
etc.)  were  used  as  initial  conditions  for  integration  in  the  corresponding  smooth 
potential.  Suppose  for  example,  that  an  initial  condition  (r0,  v0)  corresponding  to 
a chaotic  orbit  in  the  smooth  regime,  is  integrated  in  a frozen- ./V  environment  and 
that  at  time  t (where  its  coordinates  are  (r,  v)),  it  has  jumped  to  a regular  path.  If 
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Figure  4.4:  The  mean  Lyapunov  exponent  (xs),  for  ensembles  evolved  in  the 
smooth  Dehnen  potential,  selecting  as  initial  conditions  the  final  phase  space  co- 
ordinates of  orbits  that  had  been  evolved  in  frozen- iV  potentials  for  some  time  r. 

(a)  Initially  regular  (solid  line)  and  chaotic  (dashed  line)  low  energy  orbits  evolved 
with  N = 104.  (b)  Initially  regular  (solid  line)  and  chaotic  (dashed  line)  higher 

energy  orbits  evolved  with  N — 104.  (c)  The  same  as  (a)  but  for  N = 1045.  (d) 
The  same  as  (b)  but  for  N = 1045.  (e)  The  same  as  (a)  but  for  N = 105.  (f)  The 
same  as  (b)  but  for  iV  = 105. 


someone  were  to  evolve  the  initial  condition  (r,  v)  in  the  smooth  potential  he  would 
compute  a Lyapunov  exponent  that  converges  to  zero. 

The  results  of  this  experiment  are  exhibited  in  Figure  4.4  where  the  left  panels 
correspond  to  low  energies  and  the  right  panels  correspond  to  high  energies.  It 
is  obvious  that  the  average  Lyapunov  exponent  decreases  for  the  initially  chaotic 
ensembles  and  increases  for  the  initially  regular  ensembles.  The  two  lines  move 
towards  convergence.  This  convergence  happens  faster  for  smaller  N than  for  larger 
N,  which  is  a consequence  of  the  stronger  discreteness  effects  associated  with  small 
N. 
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The  fact  that  the  regular  orbits  approach  the  convergence  line  faster  is  a 
manifestation  of  the  phase  space  structure  of  the  Dehnen  potential.  Ensembles  of 
regular  orbits  correspond  to  box  orbits  which  pass  very  close  to  the  center  of  the 
system,  but  close  to  the  center  chaotic  and  regular  orbits  are  entangled  in  a very 
complex  fashion  which  makes  transitions  comparatively  easy. 


CHAPTER  5 
CHAOTIC  MIXING 


A very  interesting  question  associated  with,  but  not  limited  to,  the  question 
of  the  validity  of  the  continuum  limit  is  how  ensembles  of  initially  localised  orbits 
evolved  in  a frozen- N environment  disperse  in  the  phase  space.  This  question  was 
examined  for  three  different  types  of  ensembles:  regular,  “sticky”  chaotic,  and 
wildly  chaotic  orbits.  It  is  also  very  important  to  compare  the  behaviour  of  the 
dispersion  for  every  type  of  initial  condition  with  its  corresponding  dispersion  in 
the  smooth  potential. 

In  order  to  quantify  the  data  the  evolution  of  the  dispersion  was  computed. 

The  dispersion  is  expressed  mathematically  as: 

= (r2)  ~ (r)2  (5.1) 

where 

+ + (5.2) 

1 i=l 

The  first  potential  investigated  was  the  three-dimensional  harmonic  oscillator 
(Equation  2.3).  Since  the  smooth  potential  in  this  case  is  integrable  and  since 
for  this  particular  potential  all  orbits  oscillate  with  the  same  frequency  there  is 
no  phase  mixing  of  the  orbits  in  the  smooth  environment.  The  initial  conditions 
stay  as  localised  as  they  were  in  the  beginning  of  the  evolution.  However  in  the 
frozen- TV  analogue  of  this  potential  there  is  a dispersion  of  the  orbits  that  grows  as 
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Figure  5.1:  The  configuration  dispersion  o>  associated  with  an  initially  localised 
ensemble  evolved  in  frozen- TV  realisations  of  the  integrable  ellipsoid  potential  (2.5) 
for  variable  TV.  (a)  TV  = 1025.  (b)  TV  = 103.  (c)  TV  = 1035.  (d)  TV  = 104.  (e) 
TV  = 104-5.  (f)  TV  = 105. 


a square  root  of  time.  The  formula  that  expresses  this  evolution  is  given  as  follows: 

CTr  = (/-)*  (5.3) 

tG 

where  tG  is  the  growth  time: 

= AGNtD  (5.4) 

with  Aq  is  of  order  unity  and  to  represents  the  dynamical  time.  The  evolution  of 
a2  is  exhibited  in  Figure  5.1  for  different  values  of  TV. 

The  evolution  of  initially  localised  ensembles  of  initial  conditions  was  also 
computed  for  the  case  of  the  Plummer  potential  (Equation  2.1).  This  potential  is 
characterized  by  periodic  orbits  that  oscillate  with  different  frequencies,  so  that 
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there  is  linear  phase  mixing  even  in  the  smooth  potential  (a  oc  t).  The  time 
evolution  of  the  dispersion  was  different  for  different  values  of  N.  For  small  N the 
dispersion  grows  as  t1/2  but  for  larger  N the  dependence  becomes  linear.  This  is 
completely  understandable.  When  N is  still  small  discreteness  effects  are  quite 
strong,  so  that  a diffusion  process  (~  f1/2)  should  be  expected.  As  N increases  the 
frozen- IV  environment  approaches  the  continuum  limit  and  a behaviour  similar  to 
that  in  the  smooth  potential  (which  was  ~ t)  should  be  observed. 

The  time  evolution  of  dispersion  for  different  N can  be  seen  in  Figure  5.2.  For 
N < 104  the  square  root  law  dominates  but  for  larger  N the  behaviour  approaches 
more  and  more  the  evolution  that  characterizes  the  smooth  potential  (last  panel) 
which  is  linear.  The  spread  of  ar  is  reduced  for  smaller  N because  the  discreteness 
effects  “fuzz  out”  the  systematic  oscillations  associated  with  strictly  periodic  orbits 
in  the  unperturbed  potential. 

The  three-dimensional  harmonic  oscillator  plus  a perturbation  term  (Equation 
2.6)  is  a non-integrable  potential  where  almost  all  the  orbits  are  chaotic.  The 
dispersion  in  this  case  evolves  roughly  exponentially  in  time  and  discreteness  effects 
only  serve  to  accelerate  this  growth. 

An  example  of  the  growth  in  dispersion  for  wildly  chaotic  orbits  the  evolution 
can  be  seen  in  Figure  5.3.  The  evolution  of  the  dispersion  ar  is  not  strictly 
exponential  but  is  clearly  faster  than  the  exponential  growth  associated  with  the 
smooth  potential.  This  exponential  evolution  is  not  the  only  difference  from  the 
case  of  regular  orbits  integrated  in  the  Plummer  case.  Another  very  important 
difference  is  that  aT  for  regular  orbits  integrated  with  N as  small  as  104  5 resembles 
the  ar  evolution  in  the  corresponding  smooth  potential.  By  contrast  it  is  clear  that, 
for  the  wildly  chaotic  orbits  this  is  no  longer  true.  It  is  obvious  that  even  for  N as 
large  as  105'5  phase  mixing  happens  much  faster  than  in  the  corresponding  smooth 
potential. 
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Figure  5.2:  The  configuration  dispersion  ay  associated  with  an  initially  localised 
ensemble  evolved  in  frozen- iV  realisations  of  the  integrable  Plummer  potential  (2.2) 
for  variable  N.  (a)  N = 103.  (b)  N — 103  5.  (c)  N — 104.  (d)  N = 104  5.  (e) 
N — 105.  (f)  Evolution  in  the  smooth  potential. 
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Figure  5.3:  The  configuration  dispersion  ar  associated  with  an  initially  localised  en- 
semble of  wildly  chaotic  orbits  evolved  in  frozen- N realisations  of  the  nonintegrable 
potential  (2.6)  for  variable  N.  (a)  N — 1025.  (b)  N = 103.  (c)  N = 103  5.  (d) 
N = 104.  (e)  N = 104’5.  (f)  N = 105.  (g)  N = 105'5.  (h)  Evolution  in  the  smooth 
potential. 
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Figure  5.4:  The  configuration  dispersion  <rr  associated  with  an  initially  localised 
ensemble  of  ‘sticky’  chaotic  orbits  evolved  in  frozen- iV  realisations  of  the  noninte- 
grable  potential  (2.6)  for  variable  N.  (a)  N = 102'5.  (b)  N = 103.  (c)  N = 103'5. 
(d)  N = 104.  (e)  N = 104-5.  (f)  N = 105.  (g)  N = 105-5.  (h)  Evolution  in  the 
smooth  potential. 
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Analogous  simulations  were  also  performed  for  ensembles  of  “sticky”  chaotic 
orbits  and  the  evolution  of  their  dispersions  crr  for  different  values  of  N can  be  seen 
in  Figure  5.4.  For  N < 104  discreteness  effects  are  still  strong  and  the  “stickiness” 
is  not  clearly  manifested,  but  for  larger  N the  growth  of  a is  slower  than  for 
the  “sticky”  ensemble  exhibited  in  Figure  5.3.  It  is  important  to  note  that  the 
dispersion  saturates  in  a time  much  longer  than  the  corresponding  time  required  for 
the  wildly  chaotic  orbits. 


CHAPTER  6 

LANGEVIN  DESCRIPTION 


One  of  the  most  interesting  issues  investigated  in  this  thesis  was  to  what 
extent  it  is  possible  to  simulate  motion  in  a frozen- iV  environment  using  a Langevin 
description.  The  Langevin  description  simply  consists  of  the  Hamiltonian  equations 
of  motion  that  are  derived  for  a smooth  potential  plus  two  more  terms:  a term  that 
simulates  the  dynamical  friction  generated  physically  by  the  drag  of  the  medium, 
and  a term  that  simulates  the  random  kicks  the  test  particle  ‘feels’  because  of 
the  close  encounters  with  the  field  particles.  Therefore  the  general  form  of  the 
Langevin  equation  is: 

= -Va$  - + Fa,  a = x,y,z.  (6.1) 

where  77  can  in  principle  depend  on  r and  v but  for  the  experiments  performed 
it  had  a constant  value,  and  Fa  represents  a Gaussian  white  noise.  Fa  can  be 
completely  characterized  by  its  first  moment  and  its  correlation  function.  Its  first 
moment  is: 


(Fa(t))  = 0,  ( a,b  = x,y,z ) 

and 


{Fa(tx)Fb(t2))  = 2yQ8ab5D{tl  - t2),  (6.2) 

where  the  “temperature”  0 is  associated  with  internal  degrees  of  freedom  of  the 
system.  It  can  be  put  equal  to  the  initial  energy  of  the  initial  condition  in  order 
to  ensure  that  the  average  energy  of  the  orbit  remains  unchanged.  It  is  clear  that 
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the  noise  depends  on  the  friction  and  the  temperature  of  the  system,  which  is 
a result  of  one  of  the  most  important  theorems  in  the  statistical  mechanics,  the 
Fluctuation-Dissipation  theorem.  This  profound  theorem  connects  the  response  of 
a system  to  its  properties. 

However  since  the  energy  for  the  frozen- N orbits  remains  constant,  while 
for  the  ‘noisy’  orbits  it  does  not.  it  was  imperative  to  include  in  the  numerical 
program  an  energy  conservation  routine  that  corrected  slightly  the  orbit.  This 
routine  was  renormalising  the  velocity  at  each  time  step  by  an  overall  factor,  i.e. 
v(t  + St)  — » av(t  + St)  with  a so  chosen  that  E(t  + St)  = E(t). 

At  this  point  there  is  a hidden  technicality  associated  with  the  Langevin 
description  of  motion  in  a frozen-A’  potential.  What  happens  in  the  frozen-A" 
environments  clearly  depends  on  N.  The  Langevin  description  depends  on  the 
choice  of  77,  which  means  that  there  must  be  a logical  way  to  associate  the  right 
r]  with  the  right  N.  The  argument  we  used  to  connect  the  two  quantities  is 
straightforward.  The  value  of  7/  defines  a characteristic  relaxation  time  = I/77. 
This  is  the  time  needed  in  order  for  significant  changes  in  conserved  quantities 
like  energy  to  be  observed  (for  a self-gravitating  system).  On  the  other  hand 
it  is  well-known  that  discreteness  effects  associated  with  close  encounters  in  an 
Ar-body  system  impose  a relaxation  time  which  is:  tft  oc  ( In  N/N)tp  where  tp 
is  the  characteristic  dynamical  time.  The  last  two  formulas  give  a clear  way  to 
correspond  values  of  rj  to  values  of  N although  it  is  still  crude  and  needs  refinement 
(I  will  talk  later  about  this).  The  obvious  formula  that  connects  the  two  quantities 
is  77  oc  In  N/N  which  for  large  N becomes  77  oc  1/AL 

A good  question  before  1 proceed  is  what  is  the  law  governing  the  divergence 
between  perturbed  and  unperturbed  regular  and  chaotic  orbits  moving  in  a smooth 
potential  and  subjected  to  friction  and  white  noise.  The  behaviour  is  quite  different 
for  the  two  kinds  of  orbits:  The  divergence  of  regular  orbits  is  characterized  by 
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Figure  6.1:  The  configuration  dispersion  ar  associated  with  an  initially  localised 
ensemble  evolved  in  the  smooth  Plummer  potential  (2.2),  perturbed  by  noise  with 
0 = 1.0  and  variable  rj.  (a)  77  = 10-3.  (b)  77  = 10-3'5.  (c)  77  = 10-4.  (d)  77  = 10-4  5. 
(e)  77  = 10-5.  (f)  77  = 1CT5,5. 


a power  law.  Alternatively,  at  least  for  weak  perturbations  the  divergence  for 
chaotic  orbits  is  characterized  by  an  exponential  law  with  a rate  that  is  comparable 
to  the  largest  short-time  Lyapunov  exponent  of  the  unperturbed  orbit.  If  the 
perturbation  is  stronger  then  the  separation  between  perturbed  and  unperturbed 
orbit  quickly  becomes  macroscopic  at  which  point  the  divergence  evolves  slower 
than  exponentially. 

The  extent  to  which  a Langevin  equation  with  white  noise  can  simulate 
discreteness  effects  associated  with  motion  in  a frozen- A'  system  was  tested  in 
two  ways:  (a)  visual  comparison  of  the  plots  for  the  dispersions  aT  generated  for 
frozen- ./V  and  noisy  ensembles,  and  (b)  a comparison  of  the  slopes  associated  with 
the  growth  law  for  crr. 
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Figure  6.2:  The  configuration  dispersion  o>  associated  with  an  initially  localised  en- 
semble of  wildly  chaotic  orbits  evolved  in  the  presence  of  noise  of  the  nonintegrable 
ellipsoid  potential  (2.4)  for  0 = 1.0  and  variable  77.  (a)  77  = 10~20.  (b)  77  = 10-2'5. 
(c)  77  = 1CT30.  (d)  77  = 10~3-5.  (e)  77  = 10-40.  (f)  77  = 10_4,5.(g)  77  = 1CT50.  (h) 

77  = 10-75. 


It  is  obvious  from  Figure  6.1  (compare  it  to  Figure  5.2)  that  for  the  Plummer 
potential  there  is  a significant  cjualitative  (visual)  agreement  between  the  plots 
describing  the  time  evolution  of  ar  in  the  two  regimes,  frozen- N and  ‘noisy’,  at 
least  for  N > 104.  For  very  small  values  of  N there  are  significant  differences  since 
the  frozen- IV  ar  grows  considerably  more  rapidly  than  the  noisy  one. 

There  is  also  a considerable  similarity  for  the  case  of  chaotic  orbits,  as  is 
evident  from  Figure  6.2  (compare  it  with  Figure  5.3)  for  wildly  chaotic  orbits  and 
Figure  6.3  (compare  it  with  Figure  5.4)  for  “sticky”  chaotic  orbits.  Assuming  that 
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Figure  6.3:  The  configuration  dispersion  aT  associated  with  an  initially  localised 
ensemble  of  ‘sticky’  chaotic  orbits  evolved  in  the  presence  of  noise  of  the  nonin- 
tegrable  ellipsoid  potential  (2.4)  for  0 = 1.0  and  variable  77.  (a)  77  = 10-2  0.  (b) 
77  = 10-2'5.  (c)  77  = 10-30.  (d)  77  = 10“3'5.  (e)  77  = 10-40.  (f)  77  = 10~4'5.  (g) 
rj  = lO-50.  (h)  77  = IO-75. 


the  correspondence  between  N and  77  can  be  extrapolated  for  larger  N it  follows 
that  the  particle  number  required  for  chaotic  orbits  in  order  for  the  discreteness 
effects  to  become  unimportant  is  as  large  as  N ~ 108. 

The  exact  association  between  values  of  N and  values  of  77  can  be  accom- 
plished if  we  notice  that  the  noisy  dispersions  are  also  well  fit  by  a f1/2  growth  law 
where,  however, 


tG  = BtD/r], 


(6.3) 


with  B a constant  of  order  unity.  A comparison  of  eqs.  6.3  and  5.4  implies  a 
natural  identification: 


51 


l°gio  V = A - log10  N,  (6.4) 

with  A yet  another  constant.  For  the  ellipsoid  potential,  the  best  fit  A = O.OiO.l. 
For  the  Plummer  potential  A = 0.6±0.1.  The  identification  between  Figures  5.2 
and  6.1  was  effected  assuming  A = 0.5.  The  extent  to  which  the  growth  rates  for 
frozen- A'  and  noisy  orbits  can  be  related  by  the  simple  relation  6.4  can  be  gauged 
from  Figure  6.4,  which  superimpose  plots  of  tG(N)  and  tG{rj)  for  the  Plummer  and 
ellispoid  potentials,  with  N and  rj  related  by  eq.  6.4. 
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Figure  6.4:  The  mean  complexity  n for  orbits  evolved  in  the  potentials  (2.3), 

(2.1),  and  (2.6),  considering  both  frozen-iV  orbits  with  variable  N (diamonds)  and 
smooth  orbits  perturbed  by  noise  with  0 = 1.0  and  variable  77  (triangles).  77  was 

related  to  N by  the  relation  77  = eA/N , with  A determined  as  described  in  the  text. 
The  dashed  horizontal  line  exhibits  the  mean  complexity  of  orbits  with  the  same 
initial  conditions  evolved  in  the  smooth  potential  in  the  absence  of  noise,  (a)  Reg- 
ular orbits  in  the  Plummer  potential,  (b)  Regular  orbits  in  the  three-dimensional 
harmonic  oscillator  potential,  (c)  Wildly  chaotic  orbits  in  the  three-dimensional 
harmonic  oscillator  plus  perturbation  potential,  (d)  ‘Sticky’  chaotic  orbits  in  the 
ellipsoid  plus  black  hole  potential. 


CHAPTER  7 

NUMERICAL  ESTIMATION  OF  SMOOTH  LYAPUNOV  EXPONENTS 


There  is  one  question  left  ananswered  and  it  is  important  we  try  to  investigate 
it.  I have  mentioned  that  the  Lyapunov  exponents  computed  for  the  frozen- N 
distributions  have  values  much  larger  than  the  values  they  have  for  the  corre- 
sponding smooth  potential  integration.  The  Lyapunov  exponents  xn  associated 
with  frozen -N  descriptions  characterize  the  divergence  associated  with  microscopic 
effects,  while  the  Lyapunov  exponents  xs  computed  for  orbits  evolving  in  smooth 
potentials  are  associated  with  the  macroscopic  effect  of  the  smooth  potential.  For 
a frozen-V  orbit  the  contribution  of  the  smooth  potential  Lyapunov  exponent 
is  hidden  under  the  significantly  larger  contribution  of  the  N-body  Lyapunov 
exponent.  The  question  then  is  the  following:  is  it  possible  by  any  method  to  ex- 
tract the  smooth  Lyapunov  exponent  from  the  evolution  of  an  orbit  in  a frozen- Y 
environment? 

The  answer  to  this  question  seems  to  be  “yes”  but  only  when  N is  large 
enough  (at  least  104  5).  The  key  to  this  answer  is  given  by  the  careful  inspection 
of  the  early  evolution  of  the  divergence  between  two  orbits  starting  from  very  close 
distance  in  the  phase  space  and  evolving  in  the  same  frozen-Y  environment.  Then 
different  behaviour  is  met  depending  on  whether  the  orbit  corresponds  to  a regular 
or  a chaotic  one  in  the  smooth  environment. 

If  the  corresponding  orbit  in  the  smooth  regime  is  regular  then  the  two  orbits 
diverge  initially  exponentially  at  a rate  A ~ Xn  and  this  divergence  quickly 
saturates  once  the  separation  between  the  orbits  has  become  large  compared  to  the 
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typical  distance  between  the  point  masses.  Then  this  initial  exponential  divergence 
is  replaced  by  a power  law  divergence  which  proceeds  on  time  scale  tG  oc  Nl!2tD. 

If  the  corresponding  orbit  in  the  smooth  potential  is  chaotic  then  again  the 
initial  divergence  is  exponential  at  a rate  A ~ xn,  but  then  this  is  replaced  by  an 
exponential  divergence  at  a rate  A ~ xs,  and  eventually  it  ends  up  to  a power  law 
divergence  which  proceeds  on  time  scale  tG  oc  (In  N)tD. 

These  statements  can  become  clearer  if  one  tracks  the  evolution  of  the 


cy  . 

separation: 


5r{t)  = |rx  - r2 


(7.1) 


and  then  in  order  to  compute  the  average  divergence  for  a phase  space  region  we 
can  select  an  ensemble  of  initial  conditions  from  that  region  evolve  them  and  then 
compute  the  quantity: 

(6r)  = ~ it,  Sri  (7-2) 

Figure  7.1  exhibits  the  results  of  such  computation  for  an  ensemble  of  100 
orbits  generated  from  regular  initial  conditions  in  the  Plummer  potential.  Figure 
7.2  exhibits  the  results  for  chaotic  initial  conditions  in  the  three-dimensional  plus 
perturbation  potential.  The  N for  the  frozen- N environments  varies  from  104'5  to 
106,  and  the  initial  conditions  sampled  a phase  space  region  of  size  A r ~ Av  ~ 
10-3.  The  different  behaviour  for  regular  and  chaotic  orbits  is  obvious.  After  the 
initial  exponential  divergence  associated  with  the  microchaos  the  chaotic  orbits 
move  to  another  exponential  divergence  of  which  the  exponent  is  associated  with 
the  macrochaos,  while  the  regular  orbits  move  to  a power  law  regime. 
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Figure  7.1:  The  mean  spatial  separation  between  orbits  generated  in  the  frozen- iV 
analogue  of  the  Plummer  potential,  from  initial  conditions  seperated  in  configura- 
tion space  by  a distance  6r( 0)  = 10~5.  The  four  curves  correspond  (from  top  to 

bottom)  to  N = 104'5,  N = 105 , N — 105'5 , N = 106.  Each  curve  was  generated  by 
averaging  over  100  pairs  of  initial  conditions. 
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Figure  7.2:  The  mean  spatial  separation  between  orbits  generated  in  the  frozen- N 
analogue  of  the  three-dimensional  plus  perturbation  potential,  from  initial  condi- 
tions seperated  in  configuration  space  by  a distance  5r(0)  = 10-5.  The  four  curves 
correspond  (from  top  to  bottom)  to  N = 104'5,  N = 105,  N = 105'5,  N = 106. 
The  solid  curve  corresponds  to  a slope  of  0.022,  generated  as  a least  squares  fit. 

The  triple  dot-dashed  line  has  a slope  0.026,  equal  to  the  mean  value  of  the  smooth 
Lyapunov  exponent.  The  dashed  line  has  a slope  0.75,  equal  to  the  mean  value  of 
the  iV-body  Lyapunov  exponent  xn-  Each  curve  was  generated  by  averaging  over 
100  pairs  of  initial  conditions. 


CHAPTER  8 
CONCLUSIONS 


The  purpose  of  this  thesis  was  to  understand  the  validity  of  the  continuum 
limit  associated  with  the  smooth  potential  description  of  A-body  systems  and  to 
characterize  it  quantitatively  as  well  as  qualitatively. 

Lyapunov  exponents  associated  with  motion  of  orbits  in  frozen- A configu- 
rations have,  in  general,  values  that  are  significantly  larger  than  the  value  they 
have  for  the  same  initial  conditions  evolved  in  smooth  potential  environments. 

It  is  clear  that  the  Lyapunov  exponents  computed  numerically  for  orbits  that 
evolve  in  smooth  potentials  are  associated  with  the  form  of  the  bulk  potential. 
When  the  potential  is  integrable  these  Lyapunov  exponents  are  zero  but  when  the 
potential  is  nonintegrable  they  have  some  positive  value.  For  the  case  of  noninte- 
grable  potentials  this  positive  value  is  associated  with  the  macroscopic  chaos,  or 
macrochaos,  generated  by  the  potential  as  a whole.  The  orbits  evolving  in  frozen- A 
environments  are  characterized  by  Lyapunov  exponents  that  are  associated  mainly 
with  the  discreteness  effects,  i.e.  they  are  generated  by  a microscopic  chaos,  or 
microchaos.  This  does  not  mean  that  in  the  case  of  frozen- A'  configurations  associ- 
ated with  nonintegrable  potentials  there  is  no  contribution  to  the  destabilization  of 
the  orbit  because  of  the  action  of  the  potential  as  a whole  (macrochaos).  It  is  sim- 
ply that  this  contribution  is  small  compared  to  the  contribution  by  the  discreteness 
effects. 

Probably  the  most  striking  result  was  the  fact  that  the  Lyapunov  exponents 
of  orbits  evolved  in  frozen- A configurations  do  not  decrease  as  A increases,  even 
though  their  power  spectrum  becomes  simpler  and  approaches  the  form  it  has  when 
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the  orbit  is  evolved  in  the  smooth  potential.  This  means  that  in  the  case  of  motion 
in  an  IV-body  system  that  corresponds  to  an  integrable  potential  in  the  continuum 
limit,  for  very  large  N an  orbit  will  visually  look  completely  regular,  but  it  will 
still  be  characterized  by  a large  Lyapunov  exponent,  i.e.  it  will  be  chaotic  in  the 
sense  of  sensitivity  on  initial  conditions.  To  this  extent  Lyapunov  exponents  fail  to 
describe  the  regularity  or  chaos  in  a visual  sense,  since  x does  not  decrease  even 
when  the  ‘range’  of  chaos  gets  shorter. 

There  is  however  a theoretical  argument  that  can  explain  the  inconsistency. 
Since  in  the  frozen-iV  environments  chaos  is  generated  by  discreteness  effects,  i.e. 
the  random  interactions  between  a test  paricle  and  a collection  of  field  particles  one 
can  use  the  tidal  effects  associated  with  a pair  of  particles  separated  by  a distance 
comparable  to  the  typical  interparticle  separation  in  order  to  estimate  the  time 
necessary  for  the  growth  of  an  initial  perturbation.  The  tidal  acceleration  will  scale 
as: 

<5r  = (^r-V)a  ~ — j-Sr,  (8.1) 

where  m is  the  particle  mass  and  r the  separation.  Since 

r ~ n-1/3  - N-^Rsy,  (8.2) 

with  n the  characteristic  number  density  and  Rsys  the  size  of  the  system,  it  follows 
that  the  time  scale  associated  with  the  growth  of  an  initial  perturbation  will  be: 

U ~ 1 /y/Gp.  (8.3) 

which  is  comparable  to  the  dynamical  time  tD  and  independent  of  N.  As  N 
increases  the  individual  mass  (m  = M/N ) and  the  cube  of  the  separation  between 
particles,  n-1,  both  decrease  as  N -1  so  that  the  ratio  does  not  depend  on  N. 
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The  second  measure  of  quantification  of  chaos,  that  is  the  complexity  of 
the  power  spectrum  of  the  orbits,  agreed  completely  with  the  visual  qualitative 
description  of  the  orbit.  Orbits  evolving  in  smooth  integrable  potentials  are 
characterized  by  very  simple  power  spectra,  where  all  the  power  is  concentrated 
in  a small  number  of  frequencies.  The  evolution  in  the  corresponding  frozen- N 
configurations  gives  a more  complicated  power  spectrum  for  small  N as  well  as  a 
more  complicated  orbit  in  a visual  sense.  However  as  N increases  the  orbit  and  the 
power  spectrum  gradually  become  less  complicated  and  they  approach  the  form 
they  have  in  the  continuum  limit,  i.e.  the  smooth  case. 

For  the  nonintegrable  potentials  the  investigation  is  a little  more  difficult  be- 
cause the  power  spectrum  of  chaotic  orbits  have  a very  complex  visual  appearance 
and  a very  complex  power  spectrum  when  they  evolve  in  the  smooth  potential.  In 
the  frozen-iV  systems,  for  small  N the  power  spectrum  becomes  even  more  com- 
plex. However  as  N increases  the  visual  appearance  and  power  spectra  gradually 
approach  the  form  they  have  in  the  continuum  limit  case,  in  exactly  the  same  way 
that  regular  orbits  do. 

One  fundamental  question  to  be  answered  was  how  an  orbit  starting  from 
some  initial  condition  and  evolving  in  a smooth  potential  environment,  diverges 
from  an  orbit  starting  from  the  same  initial  condition  but  evolving  in  a frozen- TV 
analogue  of  the  smooth  potential.  For  position  and  velocity  variables  the  evolution 
of  this  divergence  was  in  general  linear  in  time  and  the  growth  time  associated  with 
this  evolution  was  oc  \/~N  for  the  regular  potentials  and  oc  IniV  for  the  chaotic  ones. 
This  means  that  regular  orbits  remain  close  to  the  smooth  characteristic  for  time 
much  longer  than  do  the  chaotic  orbits. 

In  terms  of  quantities  like  angular  momentum  that  are  invariant  in  the  smooth 
potential,  the  divergence  of  frozen- N orbits  from  the  smooth  characteristic  is  given 
as  (t/tj)1/2  with  tj  ~ y/N.  This  characterizes  a diffusive  process  which  suggests 
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that  discretensess  effects  can  be  modeled  as  white  Gaussian  noise  in  the  context  of 
a Langevin  or  Fokker-Planck  description. 

For  both  regular  and  chaotic  systems  phase  mixing  is  more  efficient  for  small 
A configurations  than  for  the  corresponding  smooth  potential  case.  As  A increases 
the  phase  mixing  evolution  resembles  more  and  more  the  phase  mixing  evolution  in 
the  smooth  potential. 

Transitions  between  regular  and  chaotic  orbits  can  and  do  arise  in  frozen- TV 
configurations,  an  effect  which  is  of  course  impossible  for  the  smooth  character- 
istics. These  transitions  are  more  common  and  happen  more  easily  for  small  A. 

As  A increases  and  approaches  the  continuum  limit  they  become  more  and  more 
difficult.  Transitions  between  sticky  and  wildly  chaotic  orbits  can  happen  in  the 
smooth  case  but  they  happen  on  a much  shorter  time  scale  for  small  A in  an 
A-body  environment. 

For  comparatively  large  A discreteness  effects  can  be  mimicked  by  energy- 
conserving  white  noise  with  amplitude  r]  oc  1/A.  This  definitely  agrees  with  the 
expectation  that  discreteness  effects  can  be  modeled  as  a sequence  of  incoherent 
binary  encounters.  The  results  of  simulations  in  frozen- A environments  and  in 
‘noisy’  configurations  agreed  in  statistical  sense.  Therefore  we  conclude  that 
noise  can  be  used  to  model  both  bulk  statistical  properties  of  orbit  ensembles 
and  qualitative  properties  of  individual  orbits.  This  means  that  the  results  of 
experiments  based  on  the  Langevin  description  (which  are  computationally  ‘cheap’) 
can  be  trusted  at  least  as  first  approximation,  to  give  important  information  about 
the  effects  of  graininess. 

However  it  is  clear  that  modeling  discretness  effects  using  the  Langevin 
description  is  not  absolutely  satisfactory.  For  small  A it  fails  to  give  similar  results 
with  the  A-body  description.  Moreover  the  A-body  Lyapunov  exponents  cannot 
be  correctly  computed  in  such  a regime.  Other  forms  of  noise,  instead  of  white 
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noise,  can  be  more  suitable  for  better  results.  This  is  a question  to  be  addressed  by 
further  research  on  the  subject. 

The  final  conclusion  of  this  thesis  is  that  at  least  in  terms  of  macroscopic 
properties  the  continuum  limit  is  valid  for  N — >-  oo.  However  convergence  towards 
this  limit  is  substantially  slower  for  density  distributions  that  correspond  in  the 
continuum  limit  to  nonintegrable  potentials  and  admit  a significant  measure  of 
chaotic  orbits,  than  for  density  distributions  that  correspond  in  the  continuum  limit 
to  integrable  potentials. 
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